The Data

knitr::include_graphics('img/landsat.png')

Load and inspect the data

landsat_file <- here('data/Landsat7.tif')

ls_1 <- raster(landsat_file)
ls_1
## class      : RasterLayer 
## band       : 1  (of  5  bands)
## dimensions : 1758, 3701, 6506358  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : -59564.57, 51465.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : /Users/logankozal/github/ESM244/labs/esm244_w2021_lab6_rasters/data/Landsat7.tif 
## names      : Landsat7 
## values     : 0, 255  (min, max)
plot(ls_1)

ls_2 <- raster(landsat_file, band = 2)
ls_3 <- raster(landsat_file, band = 3)
ls_4 <- raster(landsat_file, band = 4)

ls_stack <- raster::stack(landsat_file)
ls_stack
## class      : RasterStack 
## dimensions : 1758, 3701, 6506358, 5  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : -59564.57, 51465.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## names      : Landsat7.1, Landsat7.2, Landsat7.3, Landsat7.4, Landsat7.5 
## min values :          0,          0,          0,          0,          0 
## max values :        255,        255,        255,        255,        255

Preparing the data

# grouping cells to make it less memory intensive
ls_1 <- raster::aggregate(ls_1, fact = 3, fun = mean)
ls_2 <- raster::aggregate(ls_2, fact = 3, fun = mean)
ls_3 <- raster::aggregate(ls_3, fact = 3, fun = mean)
ls_4 <- raster::aggregate(ls_4, fact = 3, fun = mean)

ls_4
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : Landsat7 
## values     : 4.888889, 255  (min, max)
plot(ls_1, col=hcl.colors(n = 100, palette = 'Blues 2'))

plot(ls_2, col=hcl.colors(n = 100, palette = 'Greens 2'))

plot(ls_3, col=hcl.colors(n = 100, palette = 'Reds 2'))

plot(ls_4, col=hcl.colors(n = 100, palette = 'Reds 2'))

#Casey prepared this in advance

sbc_sf <- read_sf(here('data/county.shp')) %>%
  st_transform(crs(ls_1))
sbc_rast <- fasterize::fasterize(sbc_sf, ls_1, field = 'OBJECTI')
plot(sbc_rast)
writeRaster(sbc_rast, 'data/county.tif')
#mask everything not land using mask casey preprepared

sbc_rast <- raster(here('data/county.tif'))
plot(ls_3)

raster::mask(ls_3, sbc_rast) %>% plot()

ls_3 <- mask(ls_3, sbc_rast)
ls_4 <- mask(ls_4, sbc_rast)

working with rasters

Raster algebra

vec1 <- 1:5
vec1 
## [1] 1 2 3 4 5
vec1 * 2
## [1]  2  4  6  8 10
vec1 ^2
## [1]  1  4  9 16 25
ls_3 
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : Landsat7 
## values     : 8.444444, 255  (min, max)
ls_3 * 2
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : Landsat7 
## values     : 16.88889, 510  (min, max)
log(ls_3)
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 2.133509, 5.541264  (min, max)
log(ls_3); plot(log(ls_3))
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 2.133509, 5.541264  (min, max)

vec2 <- 6:10
vec1+vec2
## [1]  7  9 11 13 15
ls_3+ls_4
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 41, 498.7778  (min, max)

raster:calc()

ls_stack <- stack(ls_1, ls_2, ls_3, ls_4)
ls_mean <- raster::calc(ls_stack,fun = mean, na.rm = FALSE)

plot(ls_mean)

## Analysis # NVDI

knitr::include_graphics('img/spectrum.png')

knitr::include_graphics('img/ir_photo.jpg')

\[NDVI = \frac{NIR -RED}{NIR + RED}\]

ndvi <- (ls_4 - ls_3)/(ls_4 +ls_3) 

plot(ndvi, col=hcl.colors(100, 'Grays'))

is_forest <- function(x, thresh = .3) {
  y <- ifelse(x >= thresh, 1, NA)
  return(y)
}

forest <- calc(ndvi, fun = is_forest)
plot(forest, col = 'green')

ggplot and rasters

ndvi_df <- raster::rasterToPoints(ndvi) %>% 
  as.data.frame()
forest_df <- raster::rasterToPoints(forest) %>% 
  as.data.frame()

ggplot(data = ndvi_df, aes(x=x, y=y, fill=layer))+
  geom_raster()+
  geom_raster(data = forest_df, fill = 'green') +
  coord_sf(expand = 0)+
  scale_fill_gradient(low="black", high = 'white')+
  theme_void()+
  theme(panel.background = element_rect(fill = 'slateblue4'))
## Warning: Raster pixels are placed at uneven vertical intervals and will be
## shifted. Consider using geom_tile() instead.